# Colección de paquetes de Tidyverse
library(tidyverse)
# Estilos para ggplot2
library(ggthemes)
# Paletas de colores de RColorBrewer
library(RColorBrewer)
# Paletas de colores de viridis
library(viridisLite)
# Gráficos interactivos
library(plotly)
# Manejo de datos vectoriales
library(sf)
# Manejo de datos raster
library(terra)
# Manejo de datos raster
library(raster)
# Mapas interactivos
library(leaflet)
# Acceso a datos en GBIF
library(rgbif)
# Datos geoespaciales
library(geodata)
# Modelado de distribución de especies
library(dismo)
# Datos de elevación
library(elevatr)Nichos Ecológicos
Modelos de Nichos Ecológicos
Instalación y carga de paquetes
Datos de presencia mosquito Aedes aegypti
Se selecccionó la especie “Aedes Aegypt” dada la interesante relación de esta especie como vector primario del virus del dengue y su distribución altidudinal y latitudinal en el mundo.
# Nombre de la especie
# especie <- "Bradypus variegatus"
especie <- "Aedes aegypti"
# Consulta a GBIF
respuesta <- occ_search(
scientificName = especie,
hasCoordinate = TRUE,
hasGeospatialIssue = FALSE,
limit = 10000,
curlopts = list(timeout = 120)
)
# Extraer datos de presencia
presencia <- respuesta$data# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')# Leer los datos de presencia de un archivo CSV
#| label: Datos-Pharomachrus
#| warning: false
#| message: false
presencia <- read_csv('presencia.csv')Rows: 10000 Columns: 175
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (130): scientificName, issues, datasetKey, publishingOrgKey, installati...
dbl (28): key, decimalLatitude, decimalLongitude, crawlId, taxonKey, kingd...
lgl (11): isSequenced, isInCluster, organismName, gadm.level0.gid, gadm.le...
dttm (5): lastCrawled, lastParsed, dateIdentified, modified, lastInterpreted
date (1): georeferencedDate
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.presencia <- st_as_sf(
presencia,
coords = c("decimalLongitude", "decimalLatitude"),
remove = FALSE, # conservar las columnas de las coordenadas
crs = 4326
)Gráfico de barras registro por país
De acuerdo con la Organización Mundial de la Salud (OMS), el mosquito Aedes aegypti se encuentra en regiones tropicales y subtropicales de Asia Sudoriental, el Pacífico Occidental, el Mediterráneo Oriental, África y las Américas. En América, se extiende desde el sur de los Estados Unidos hasta el centro de Argentina y Uruguay.OMS (2024) Los resultado del gráfico de barras muestran a Estados Unidos como el país con mayor cantidad de registros lo cual puede estar principalmente vinculado a su infraestrutura de monitoreo, investigación científica y recursos en general.
# Filtrar los datos para incluir solo los N países con más registros
top_n <- 20 # Cambiar este valor según el número de países deseados
datos_filtrados <- presencia |>
count(countryCode, sort = TRUE) |>
head(top_n)
# Ordenar las barras de manera descendente
datos_filtrados <- datos_filtrados |>
mutate(countryCode = fct_reorder(countryCode, n, .desc = TRUE))
# Gráfico ggplot2 con datos filtrados
grafico_ggplot2 <-
datos_filtrados |>
ggplot(aes(x = fct_infreq(countryCode), y = n)) +
geom_bar(
stat = "identity",
aes(
text = paste0(
"Cantidad de registros de presencia: ", n
)
)
) +
ggtitle("Cantidad de registros de presencia por país") +
xlab("País") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: GBIF") +
theme_economist()
# Convertir el gráfico a plotly para interactividad
grafico_interactivo <- ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')
# Mostrar el gráfico interactivo
grafico_interactivoGráfico de registro por año
Según los resultados del gráfico de registros por año entre el 2017 y el 2024 existe un descenso constante con la excepción del 2020 en donde se vio un aumento de registros llegando a 2483 observaciones. Segun la OMS el creciente aumento en casos de dengue tiene una relación directa con cambios en la distribución de los vectores (principalmente los mosquitos Aedes aegypti y Aedes albopictus), especialmente en los países que no habían padecido dengue; escomo como consecuencias del fenómeno de El Niño en 2023 y el cambio climático, que incrementa las temperaturas, las precipitaciones y la humedad. OMS (2024)
# Gráfico ggplot2
grafico_ggplot2 <-
presencia |>
st_drop_geometry() |>
group_by(year) |>
summarize(n = n()) |>
ggplot(aes(x = year, y = n)) +
geom_line() +
geom_point(
aes(
text = paste0(
"Año: ", year, "\n",
"Cantidad de registros: ", n
)
)
) +
ggtitle("Cantidad de registros de presencia por año") +
xlab("Año") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: GBIF") +
theme_economist()
# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')Mapa de Registros de Aedes aegypti
Como se mencionó anteriormente el mosquito Aedes aegypti se encuentra en regiones tropicales y subtropicales de Asia Sudoriental, el Pacífico Occidental, el Mediterráneo Oriental, África y las Américas.El mapa de registros confirma estos datos y se puede notar focos principalment en México y Centro America, costa Atlántica de Brasil y el sudeste Asiático.
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Aedes aegypti"
) |>
addLayersControl(
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c("Registros de Aedes aegypti"))Carga de datos ambientales
Datos de WorldClim
# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())
# Nombres de las variables climáticas
names(clima) [1] "wc2.1_10m_bio_1" "wc2.1_10m_bio_2" "wc2.1_10m_bio_3" "wc2.1_10m_bio_4"
[5] "wc2.1_10m_bio_5" "wc2.1_10m_bio_6" "wc2.1_10m_bio_7" "wc2.1_10m_bio_8"
[9] "wc2.1_10m_bio_9" "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
[13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
[17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"Datos de elevación de Shuttle Radar Topography Mission (SRTM)
# Consulta a SRTM
elevacion <- elevation_global(res = 10, path = tempdir())
# Nombres de las variables de elevación
names(elevacion)[1] "wc2.1_10m_elev"Definición del área de estudio
# Definir la extensión del área de estudio
area_estudio <- ext(
min(presencia$decimalLongitude) - 5,
max(presencia$decimalLongitude) + 5,
min(presencia$decimalLatitude) - 5,
max(presencia$decimalLatitude) + 5
)
# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)
elevacion <- crop(elevacion, area_estudio)Mapa de registros combinado con mapas de variables ambientales
Se seleccionaron como variables ambientales la precipitación, temperatura y elevación como factores incidentes en la idoneidad del hábitat para esta especie de mosquito.
# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
palette = rev(brewer.pal(11, "RdYlBu")),
values(clima$wc2.1_10m_bio_1),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
palette = "Blues",
values(clima$wc2.1_10m_bio_12),
na.color = "transparent"
)
# Paleta de colores de elevación
colores_elevacion <- colorNumeric(
palette = "Greens", # Cambia según tu preferencia
values(elevacion), # Usa el objeto 'elevacion'
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
clima$wc2.1_10m_bio_1,
colors = colores_temperatura,
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
clima$wc2.1_10m_bio_12,
colors = colores_precipitacion,
opacity = 0.6,
group = "Precipitación",
) |>
addRasterImage( # capa raster de elevación
elevacion,
colors = colores_elevacion, # Usa el objeto 'elevacion'
opacity = 0.6,
group = "Elevación",
) |>
addCircleMarkers(
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Aedes aegypti"
) |>
addLegend(
title = "Temperatura",
values = values(clima$wc2.1_10m_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(clima$wc2.1_10m_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Elevación",
values = values(elevacion), # Usa el objeto 'elevacion'
pal = colores_elevacion,
position = "bottomleft",
group = "Elevación"
) |>
addLayersControl(
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c("Temperatura", "Precipitación", "Elevación", "Registros de Bradypus variegatus")
) |>
hideGroup("Precipitación") |>
hideGroup("Elevación")Modelización
El proceso de generación del modelo comienza con la preparación de los datos de presencia de la especie, para lo cual se extraen las coordenadas de longitud y latitud y se almacenan en un data.frame. Posteriormente, se eliminan las coordenadas duplicadas para garantizar que no haya sesgos en el modelo debido a datos redundantes. Posterior a ello, se realiza una división aleatoria de los datos en dos subconjuntos: un 70% de los datos se utiliza para entrenar el modelo, mientras que el 30% restante se reserva para evaluar su desempeño. Esta selección aleatoria se realiza de forma reproducible estableciendo la función set.seed().
Las variables climáticas necesarias para el modelado se convierten al formato ráster compatible con el paquete dismo, agrupándolas en una pila de capas (stack). Una vez preparados los datos, se entrena el modelo MaxEnt (Maximum Entropy), utilizando las coordenadas de presencia y las capas climáticas. Este algoritmo identifica patrones en las variables ambientales para estimar la idoneidad del hábitat de la especie en función de las condiciones climáticas observadas en los registros de presencia. Dicho modelo es una herramienta ampliamente utilizada en la modelización de nichos ecológicos y la predicción de distribuciones de especies (Mafuwe et al., 2022). Este enfoque aplica técnicas de aprendizaje automático para estimar la distribución geográfica potencial de una especie basándose únicamente en datos de presencia y variables ambientales, donde el principio subyacente al asumir la mínima información posible (máxima entropía), se puede modelar la distribución más uniforme que se ajuste a las restricciones impuestas por los datos disponibles (Elith et al., 2011).
Por último, el modelo entrenado se aplica a las variables climáticas para generar el mapa continuo de idoneidad del hábitat. Este mapa indica las áreas más propensas a ser adecuadas para la especie según el escenario climático, donde valores altos representan una mayor probabilidad de idoneidad. Este proceso no solo permite identificar las áreas clave para la conservación, sino que también establece una base para evaluar la distribución potencial de la especie bajo diferentes escenarios ambientales o climáticos futuros.
Creación de conjuntos de entrenamiento y de evaluación
# Crear dataframe con columnas de longitud y latitud
coordenadas_presencia <- data.frame(
decimalLongitude = presencia$decimalLongitude,
decimalLatitude = presencia$decimalLatitude
)
# Eliminar coordenadas duplicadas
coordenadas_presencia <- unique(coordenadas_presencia)
# Establecer una "semilla" para garantizar que la selección aleatoria sea reproducible
set.seed(123)
# Cantidad de registros de presencia
n_presencia <- nrow(coordenadas_presencia)
# Con sample(), se selecciona aleatoriamente una proporción (ej. 0.7)
# de los índices de los datos de presencia para el conjunto de entrenamiento
indices_entrenamiento <- sample(
1:n_presencia,
size = round(0.7 * n_presencia)
)
# Crear el subconjunto de entrenamiento utilizando los índices seleccionados
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
# Crear el subconjunto de evaluación con los datos restantes
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]
# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima <- raster::stack(clima)
# Ejecutar el modelo
modelo_maxent <- maxent(x = clima, p = entrenamiento)
# Aplicar el modelo entrenado a las variables climáticas
# para generar un mapa de idoneidad del hábitat
prediccion <- predict(modelo_maxent, clima)Evaluación
# terra::extract() extrae los valores del raster de predicción
# en las coordenadas de evaluación
# eval_pres almacena los valores de idoneidad predichos
# en los puntos de evaluación de presencia
eval_pres <- terra::extract(
prediccion,
evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)
# Generar puntos aleatorios dentro del área de estudio definida.
# Estos puntos se asumen como ausencias de la especie.
ausencias <- randomPoints(mask = clima, n = 1000)
# eval_aus almacena los valores de idoneidad predichos
# en los puntos de ausencia
eval_aus <- terra::extract(
prediccion,
ausencias
)
# Generar estadísticas de evaluación del modelo
resultado_evaluacion <- evaluate(p = eval_pres, a = eval_aus)Generación de Curva ROC
La curva ROC presenta la evaluación del desempeño del modelo binario en términos de su capacidad para diferenciar entre clases positivas y negativas, en este caso para le modelo MaxEnt, antes descrito. El área bajo la curva (AUC) es de 0.927, lo que indica un rendimiento muy satisfactorio, ya que es un valor muy cercano a 1, lo cual representa un modelo ideal. A su vez, dicho valor de AUC sugiere que el modelo tiene una gran capacidad para discriminar entre las clases.
La forma de la curva ROC, está muy por encima de la línea diagonal roja que representa un modelo aleatorio (AUC = 0.5), confirmando que el modelo es significativamente mejor que un enfoque de clasificación aleatorio. La curva se acerca al vértice superior izquierdo del gráfico, lo que indica que el modelo logra una alta tasa de verdaderos positivos (TPR) mientras mantiene una baja tasa de falsos positivos (FPR). En términos prácticos, esto significa que el modelo identifica correctamente la mayoría de los casos positivos mientras minimiza los errores al clasificar casos negativos como positivos.
Por ende, a nivel general, la forma de la curva y el alto valor del AUC demuestran que el modelo es bastante efectivo y confiable para distinguir entre las dos clases, lo que lo convierte en una herramienta adecuada para obtener una alta precisión en la clasificación de la distribución de la especie en cuestión (Aedes aegypti) mediante el modelo mencionado.
# Datos para graficar la curva ROC
datos_roc <- data.frame(
FPR = resultado_evaluacion@FPR,
TPR = resultado_evaluacion@TPR,
Umbral = resultado_evaluacion@t
)
# Valor AUC
auc <- resultado_evaluacion@auc
# Gráfico ggplot2
grafico_ggplot2 <-
ggplot(
datos_roc,
aes(
x = FPR,
y = TPR,
u = Umbral
)
) +
geom_line(
color = "blue",
size = 1
) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
labs(title = paste("Curva ROC (AUC =", round(auc, 3), ")"),
x = "Tasa de falsos positivos (FPR)",
y = "Tasa de verdaderos positivos (TPR)") +
theme_minimal()
# Gráfico plotly
ggplotly(grafico_ggplot2) |>
config(locale = 'es')Mapa continuo
Se utilizaron como variables climáticas la precipitación, temperatura y elevación para predecir a partir del modelo la idoneidad del hábitat para esta especie. Como resultado del modelo se puede notar una distribución más amplia alrededor de la ubicación de los registros, extendiéndose más al norte de los Estados Unidos, hacia el interior de Brasil, Paraguay y Bolivia, al norte del sudeste asiático, las costas africanas y al este de Australia.
Si bien es cierto el modelo se utilizó para definir aspecto de idoneidad de hábitat, no pareció interesante como identificación de estos espacios pueden incidir en temas de salud humana, en este caso como vector de enfermedades.
# Paleta de colores del modelo
colores_modelo <- colorNumeric(
palette = c("white", "black"),
values(prediccion),
na.color = "transparent"
)
# Crear la paleta de colores para la elevación
colores_elevacion <- colorNumeric(
palette = terrain.colors(10), # Paleta para elevación
domain = values(elevacion),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
clima$wc2.1_10m_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura"
) |>
addRasterImage( # capa raster de precipitación
clima$wc2.1_10m_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación"
) |>
addRasterImage( # capa raster de elevación
elevacion,
colors = colores_elevacion,
opacity = 0.6,
group = "Elevación"
) |>
addRasterImage( # capa raster del modelo de distribución
prediccion,
colors = colores_modelo,
opacity = 0.6,
group = "Modelo de distribución"
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = presencia,
stroke = F,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Aedes aegypti"
) |>
addLegend(
title = "Temperatura",
values = values(clima$wc2.1_10m_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(clima$wc2.1_10m_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Elevación (m)",
values = values(elevacion),
pal = colores_elevacion,
position = "bottomleft",
group = "Elevación"
) |>
addLegend(
title = "Modelo de distribución",
values = values(prediccion),
pal = colores_modelo,
position = "bottomright",
group = "Modelo de distribución"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Temperatura",
"Precipitación",
"Elevación",
"Modelo de distribución",
"Registros de Aedes aegypti"
)
) |>
hideGroup("Temperatura") |>
hideGroup("Precipitación") |>
hideGroup("Elevación")Mapa binario
El resultado del mapa vinario permite apreciar puntual la distribución espacial de la especie, sin que se consideren la cantidad de observaciones, sino más bien la presencia o ausencia de las mismas. Se continua con focos en la costa de Brasil, Centro América y México, el sudeste Asiático y a pesar que Australia únicamente refleja un cerca de tres observaciones para este análisis se vuelve mucho más visible.
# Definir el umbral
umbral <- 0.5
# Crear el raster binario
prediccion_binaria <- (prediccion >= umbral) * 1
# Crear la paleta de colores para el raster binario
colores_prediccion_binaria <- colorFactor(
palette = c("transparent", "blue"), # "transparent" para las áreas no adecuadas
domain = c(0, 1),
na.color = "transparent"
)
# Crear la paleta de colores para elevación
colores_elevacion <- colorNumeric(
palette = terrain.colors(10), # Paleta para elevación
domain = values(elevacion),
na.color = "transparent"
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage(
prediccion_binaria,
colors = colores_prediccion_binaria,
opacity = 0.6,
group = "Modelo de distribución binario"
) |>
addRasterImage(
elevacion,
colors = colores_elevacion,
opacity = 0.6,
group = "Elevación"
) |>
addCircleMarkers(
data = presencia,
stroke = FALSE,
radius = 3,
fillColor = 'red',
fillOpacity = 1,
popup = paste(
paste0("<strong>País: </strong>", presencia$country),
paste0("<strong>Localidad: </strong>", presencia$locality),
paste0("<strong>Fecha: </strong>", presencia$eventDate),
paste0("<strong>Fuente: </strong>", presencia$institutionCode),
paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de Bradypus variegatus"
) |>
addLegend(
title = "Modelo de distribución binario",
labels = c("Ausencia", "Presencia"),
colors = c("transparent", "blue"),
position = "bottomright",
group = "Modelo de distribución binario"
) |>
addLegend(
title = "Elevación (m)",
values = values(elevacion),
pal = colores_elevacion,
position = "bottomright",
group = "Elevación"
) |>
addLayersControl(
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Modelo de distribución binario",
"Elevación",
"Registros de Aedes aegypti"
)
) |>
hideGroup("Elevación")